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Abstract 


We  present  a  summary  of  investigations  on  the  use  of  proper  orthogonal  decomposition 
(POD)  techniques  as  a  reduced  basis  method  for  computation  of  feedback  controls  and 
compensators  in  a  high  pressure  chemical  vapor  deposition  (HPCVD)  reactor  that  in¬ 
cludes  multiple  species  and  controls,  gas  phase  reactions,  and  time  dependent  tracking 
signals  that  are  consistent  with  pulsed  vapor  reactant  inputs.  Numerical  implementa¬ 
tion  of  the  model-based  feedback  control  uses  a  reduced  order  state  estimator,  based  on 
partial  state  observations  of  the  fluxes  of  reactants  at  the  substrate  center,  which  can 
be  achieved  with  current  sensing  technology.  We  demonstrate  that  the  reduced  order 
state  estimator  or  compensator  system  is  capable  of  substantial  control  authority  when 
applied  to  the  full  system. 
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1  Introduction 


Stringent  control  of  layer  thickness  and  composition  is  essential  in  the  production  of 
advanced  optoelectronic  integrated  circuits.  This  can  sometimes  be  addressed,  in  part, 
through  open-loop  optimization  [7,  25,  41].  However,  process  variability  and  the  in¬ 
creasing  demands  put  on  control  of  layer  thickness  and  composition  in  state-of-the  art 
devices  increase  the  desirability  of  real-time  control  of  film  growth  [8,  19,  21,  43,  42], 
Low  pressure  chemical  vapor  deposition  (CVD)  processes  are  the  preferred  choice  for 
manufacturing  many  advanced  optoelectronic  devices.  On  the  other  hand,  there  are  also 
materials  of  potential  industrial  use  ( e.g .,  InN  films)  that  exhibit  relatively  high  decom¬ 
position  pressure  as  compared  to  other  III-V  compounds.  These  can  not  be  adequately 
produced  at  desirable  process  temperatures  under  low  pressure  conditions  and  thus  it  is 
be  desirable  to  extend  the  CVD  processing  to  higher  pressures.  We  are  presently  col¬ 
laborating  with  material  scientists  at  N.C.  State  University  to  design  and  build  such  a 
HPCVD  reactor  with  real-time  sensing  and  control  as  an  innovative  feature  of  this  proto¬ 
type  reactor.  Previous  work  within  this  collaboration  has  experimentally  demonstrated 
the  successful  implementation  of  closed-loop  control  of  film  thickness  and  composition  in 
GaP/Gai_xInxP  heterostructures  grown  in  a  low  pressure  pulsed  chemical  beam  epitaxy 
(CBE)  reactor  utilizing  real-time  optical  p-polarized  reflectance  signal  (PRS)  sensing  [19]. 

The  modeling  of  reactant  transport  in  HPCVD  reactors  is  more  complicated  than  in 
the  case  of  low  pressure  CVD  reactors.  In  general,  a  full  mathematical  model  describ¬ 
ing  transport  processes  in  (HPCVD)  systems  is  given  by  a  system  of  nonlinear  partial 
differential  equations  representing  the  continuity,  momentum,  energy,  and  species  equa¬ 
tions  of  state.  Numerical  simulations  and  control  designs  of  such  systems  using  finite 
element,  finite  difference,  or  spectral  methods  lead  to  very  large  systems  of  ordinary 
differential  equations  rendering  real-time  full  model-based  feedback  control  design  infea¬ 
sible.  Substantial  dimensional  reduction  can  be  realized  using  the  method  of  POD  to 
more  efficiently  represent  the  system  data  [33].  In  recent  work,  we  demonstrated  the 
potential  of  a  POD  based  reduced  order  model  as  a  basis  for  implementing  real  time 
control  by  presenting  a  method  for  the  design  of  state  feedback  controllers  based  on  data 
provided  by  accurate  model  simulations  of  species  transport  [26] . 

In  this  paper  we  extend  the  design  of  state  feedback  controllers  to  include  experimentally 
relevant  conditions  of  high  pressure  III-V  film  growth  with  multiple  species  and  controls, 
gas  phase  reactions,  and  time  dependent  tracking  signals  that  are  consistent  with  pulsed 
vapor  reactant  inputs.  We  present  implementation  of  Dirichlet  boundary  control  of  dilute 
Group  III  and  Group  V  reactants  transported  by  convection  and  diffusion  to  an  absorbing 
substrate  while  undergoing  gas  phase  reactions.  Computational  fluid  dynamics  (CFD) 
simulations  provide  data  on  the  full  system  and  are  used  to  construct  a  POD  reduced 
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order  model. 


The  control/compensator  problem  (or  LQG  tracking  problem)  is  formulated  as  a  linear 
quadratic  regulator  (LQR)  time  dependent  tracking  problem.  The  state  of  the  system  is 
estimated  via  a  compensator  gain  from  observations  of  the  fluxes  of  Group  III  and  Group 
V  reactants  to  the  substrate,  which  are  assumed  to  be  proportional  to  the  growth  rate, 
with  the  observed  fluxes  attempting  to  track  time  dependent  target  flux  values  such  as 
is  encountered  in  pulsed  CVD  him  growth.  For  example,  one  might  choose  a  target  flux 
such  that  the  integrated  pulse  deposits  enough  material  to  form  a  monolayer.  Dirichlet 
boundary  control  at  the  inlet  determines  the  mass  fractions  of  the  incoming  Group  III 
and  Group  V  reactants.  We  demonstrate  that  the  reduced  order  state  estimator  or 
compensator  system  is  capable  of  substantial  control  authority  when  applied  to  the  full 
system. 

The  use  of  POD  (also  known  under  the  names  of  principal  component  analysis  [23]  and 
Karhunen-Loeve  expansion  [24,  30])  as  a  method  of  feature  extraction  from  data  sets  is 
well  known  in  statistical  and  pattern  recognition  fields  [20]  and  has  been  applied  in  a  wide 
variety  of  areas,  such  as  materials  processing  [33,  41],  characterization  of  human  faces 
[37],  and  turbulent  coherent  Hows  ([4,  15,  16,  18,  22,  31,  38]  -  see  also  the  surveys  [14,  32]). 
The  POD  method  is  a  linear  transformation  of  a  multivariate  data  set  into  an  optimal  set 
of  uncorrelated  variables  (POD  modes).  The  original  multivariate  data  can  be  written  as 
linear  combinations  of  the  POD  modes.  In  many  cases  the  POD  modes  more  efficiently 
describe  the  variability  of  the  original  data  and  some  dimensional  reduction  is  possible 
by  retaining  only  the  most  important  modes. 

Recent  use  of  POD  for  reduction  of  order  in  distributed  parameter  systems  includes 
applications  to  parameter  estimation  or  inverse  problems  [11]  and  both  open  loop  and 
feedback  control  design  [1,  2,  9,  10,  26,  27,  33,  34,  41],  There  are  a  number  of  nontrivial 
issues  related  to  the  use  of  reduced  basis  methods  in  general,  and  POD  methods  in  par¬ 
ticular,  as  a  foundation  for  approximation  methods  in  infinite  dimensional  systems  such 
as  those  modeled  by  distributed  parameter  or  partial  differential  equation  systems.  The 
most  important  issue  is  whether  the  infinite  dimensional  system  can  be  approximated 
well  by  a  finite  span  of  basis  elements.  Initial  efforts  on  error  analysis  for  POD  based 
approximations  in  forward  or  simulation  problems  can  be  found  in  [28].  A  more  fun¬ 
damental  question  in  regard  to  order  reduction  is  whether  the  important  features  of  the 
system  are  essentially  low  finite  dimensional  in  nature?  There  is  mounting  computational 
evidence  that  the  answer  to  this  question  is  positive  for  many  structural  systems  and  for 
a  substantial  number  of  fluid  and  electromagnetic  applications  [11].  However,  even  when 
this  question  can  be  answered  in  the  affirmative,  it  is  not  at  all  clear  that  one  can  use 
such  an  approach  for  control  design. 
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Specifically,  if  one  uses  ‘snapshots’  of  the  uncontrolled  system  (as  was  done  in  [9,  10]) 
to  construct  the  POD  basis  elements,  there  is  little  reason  to  expect  that  control  design 
based  on  this  finite  dimensional  approximation  will  be  effective  when  applied  to  the 
original  system.  The  controlled  original  system  itself  may  require  a  different  set  of  finite 
dimensional  elements  for  efficient  approximation  or,  worse  yet,  may  not  be  amenable  to 
a  low  order  basis  approximation.  Fortunately,  the  results  of  [9,  10]  suggest  that,  at  least 
in  some  situations,  this  is  not  the  case.  Another  approach  is  to  ‘snapshot’  on  the  system 
under  several  levels  of  nontrivial  control  inputs  that  are  not,  in  general,  derived  from  any 
optimal  or  suboptimal  design.  This  approach  has  proved  effective  in  open  loop  [33]  as 
well  as  closed  loop  [26]  control  problems  and  is  employed  in  this  paper. 

In  Section  2  below  we  describe  the  particular  system  under  investigation  here  and  describe 
the  simulations  needed  to  develop  POD  based  control  design.  A  sample  of  some  of  our 
computational  results  based  on  this  approach  is  given  in  Section  3.  Brief  conclusions  are 
then  summarized. 


2  Methods 

We  consider  a  2D  rectangular  domain  (Fig.  1)  representing  the  longitudinal  cross  section 
through  the  center  of  an  HPCVD  reactor  with  the  following  dimensions:  0.011  m  height, 
0.156  m  total  length,  and  0.048  m  substrate  length.  Gas  flow  enters  from  the  inlet  on 
the  left  of  the  domain  (Pi)  and  exits  at  the  right  boundary  of  the  domain  (r3),  after 
passing  over  the  heated  substrate  on  the  bottom  wall  (r.5).  We  control  the  concentration 
of  the  source  material  at  the  inlet  in  an  effort  to  obtain  a  desired  flux  of  reactants  to  the 
substrate,  which  can  be  assumed  to  be  proportional  to  the  growth  rate.  The  remainder  of 


Figure  1:  Two-dimensional  cross  section  of  HPCVD  reactor  geometry  with  the  following  dimensions 
height=0.011  m,  length=0.156  m,  substrate  length=0.048  m. 

the  Methods  section  is  organized  into  two  parts.  The  first  presents  methods  for  obtaining 
simulated  experimental  data  or  time  snapshots  necessary  for  obtaining  the  POD  modes 
for  the  reduced  order  model,  while  the  second  section  presents  the  control  equations. 
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2.1  Simulated  data 


We  consider  only  trace  amounts  of  reactants  mixed  with  the  carrier  gas.  Under  this 
dilute  approximation,  we  can  obtain  steady-state  solutions  for  the  velocity,  temperature, 
and  pressure  that  are  independent  of  the  reactant  concentration  and  are  governed  by  the 
following  set  of  equations 
(continuity) 


V  •  (pv)  =  0, 


(1) 


(momentum) 


pv  •  W  =  —  VP  +  V  •  r  —  pg, 


(2) 


where  the  viscous  stress  tensor  is  of  the  form 


t  =  — -/i(V  •  v)I  +  p( W  +  Vvr), 
o 


(3) 


(energy) 


pcpv  •  VT  =  V  •  (k'VT), 


(4) 


where  g  is  the  gravitational  acceleration,  u,  T,  and  P  are  the  velocity,  temperature,  and 
pressure,  p,  cp,  and  k  are  the  viscosity,  specific  heat,  and  thermal  conductivity  of  the 
carrier  gas.  The  density  variations  are  modeled  as  p  =  po[l  —  (3 (T  —  To)],  where  T0  is  a 
reference  temperature,  po  is  a  reference  density  calculated  from  the  ideal  gas  law  at  the 
reference  temperature  and  reactor  pressure,  and  [3  is  the  volume  coefficient  of  expansion 
(f 3  =  1/T).  The  boundary  conditions  for  the  above  system  of  equations  (l)-(4)  will  be 
discussed  in  subsequent  sections. 


The  governing  equations  are  discretized  using  the  Galerkin  finite  element  method  with 
weighted  residuals  for  the  degrees  of  freedom  ( v. ,  P,  and  T ).  A  mixed  formulation  with 
132  quadrilateral  elements  (corresponding  to  453  nodes)  is  used  with  piecewise  linear  dis¬ 
continuous  elements  for  pressure,  and  quadratic  (8-noded)  elements  for  the  other  degrees 
of  freedom.  Solutions  for  the  variables  v.  P,  T,  and  Yn  are  obtained  using  simulations 
with  commercially  available  code  (FIDAP,  Fluid  Dynamics  International,  Evanston,  IL) 
on  a  Silicon  Graphics  Power  Indigo  2. 


For  the  results  presented  in  this  article,  we  consider  a  hydrogen  carrier  gas  at  atmospheric 
pressure.  Temperature  dependent  values  for  p,  k,  and  cp  are  linearly  interpolated  from 
measurements  taken  from  the  available  literature  [29,  39,  13].  A  parabolic  velocity  flow 
profile  is  specified  at  the  inlet  (Tl),  with  an  average  inlet  velocity  of  0.1147  m/s.  No 
slip  (zero  velocity)  boundary  conditions  are  imposed  on  those  portions  of  the  model 
corresponding  to  the  reactor  walls  (r2,  T4,  T5,  and  T6).  Room  temperature  boundary 
conditions  are  imposed  at  the  inlet  and  along  the  upper  wall  (T2).  Along  the  bottom 
wall,  the  substrate  (r5)  temperature  is  fixed  at  800°K  ,  with  a  non-linear  temperature 
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decrease  from  the  substrate  edge  to  the  inlet  (T6)  and,  similarly,  from  the  substrate  edge 
to  the  outlet  (T4)  (see  Fig.  1). 

Steady-state  solutions  for  v.  T,  and  p  obtained  above  are  used  in  the  time-dependent 
governing  equation  for  species  transport 

+  v  •  vvn  =  -  v  •  [PDn vvn)  +  jy-m ,  (5) 

where  Dn  is  the  diffusivity  of  the  species,  Yn  is  the  mass  fraction  of  the  nth  species,  NR 
is  the  number  of  gas  phase  reactions,  and  rni  is  the  rate  of  production  of  species  n  in 
the  itli  chemical  reaction.  Solutions  of  this  equation  (with  boundary  conditions  to  be 
specified)  are  used  to  generate  data  for  construction  of  the  POD  modes. 

Trimethylindium  (TMI)  and  phosphine  are  source  materials  commonly  used  for  CVD 
growth  of  III-V  films  where  indium  and  phosphorus  are  the  constituents.  For  high  pres¬ 
sure  applications,  the  source  gases  are  introduced  in  separate  pulses  with  sufficient  time 
between  pulses  to  allow  the  reactants  to  exit  before  introduction  of  the  next  source  mate¬ 
rial  (Fig.  2).  Pulsing  of  the  III-V  source  materials  prevents  nucleation  of  the  film  in  the 
gas  phase  and  makes  PRS  observation  and  analysis  possible.  A  constant  total  volumetric 
flow  rate  is  maintained  throughout  the  pulsing  cycle. 


Figure  2:  Schematic  representation  of  source  gas  pulsing. 


Under  the  reactor  conditions  considered  here  (H2  carrier  gas,  800°K  substrate  temper¬ 
ature,  and  1  atm  pressure),  there  are  no  effective  gas  phase  reaction  mechanisms  for 
phosphine,  and  the  only  significant  gas  phase  reaction  for  TMI  is  the  decomposition  of 
TMI  to  monomethylindium  (MMI)  and  two  methyl  molecules,  In(CH3)3  — >  I11CH3  + 
2CH3.  This  reaction  can  be  described  as  a  first-order  Arrhenius  reaction 

rn  =  ’  V 7  a// •  (6) 

Wtmi 

where  vn  refers  to  the  stoichiometry  of  species  n  in  the  reaction,  Wn  and  Wtmi  refer  to 
the  molecular  weight  of  species  n  and  TMI,  respectively,  ko  =  5.25  x  1015  s-1  is  the  rate 
constant,  and  E  =  47.2  kcal/mol  is  the  activation  energy  [40]. 
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A  transient  simulation  of  TMI  transport  (5)  through  one  cycle  of  pulsing  and  clearing  - 
including  the  thermal  decomposition  of  TMI,  diffusion  of  TMI  and  MMI,  and  absorption 
at  the  substrate  -  is  used  to  obtain  data  for  construction  of  the  POD  modes  for  TMI 
and  MMI.  Similarly,  simulation  data  of  phosphine  transport  from  one  pulsing  and  clear¬ 
ing  cycle  is  used  to  construct  the  POD  modes  for  phosphine.  The  simulation  system, 
including  a  nonzero  boundary  input,  is  given  by 


+  v-  XYn 
Yn(  0,x) 
Y1(t,x),Y2(t,x) 
Y3(t,x ) 
Yn(t,  x) 
dYn(t,x) 
Bn 

n 


iv  •  (pDnXYn)  +  \n(T)Yi 
0 

1  or  0  (pulsed)  on  T1 
0  on  n 
0  on  T5 

0  on  T2  U  T3  U  T4  U  T6 
1,2,3 


(7) 


where  Yi,Y2,  and  Y3  refer  to  the  mass  fractions  of  TMI,  phosphine,  and  MMI,  respectively, 
A i(T)  =  -k0e(-E/RT\  A 2(T)  =  0,  A 3(T)  =  (Wa/W^koe^/^,  and  ^  denotes  the 
outward  normal  derivative.  Since  the  methyl  (CH3)  molecules  do  not  participate  in  film 
growth  or  otherwise  affect  the  transport  properties  (under  the  dilute  approximation), 
we  do  not  include  them  in  the  transport  equations.  The  reactor  walls  (T2,  T4,  and  T6) 
are  assumed  non- absorbing,  and  the  substrate  (r5)  is  assumed  to  be  perfectly  absorbing 
(concentration  of  zero).  Values  for  v  and  p  appearing  in  (5)  are  provided  by  the  steady- 
state  solutions.  Temperature  dependent  values  for  the  diffusivities  Dn  in  hydrogen  are 
linearly  interpolated  from  values  taken  from  the  available  literature  [40]. 


Initially,  the  reactant  mass  fractions  are  assumed  zero  everywhere  in  the  reactor  and  only 
the  H2  carrier  gas  is  flowing.  At  the  start  of  the  pulsing  cycle,  a  normalized  TMI  mass 
fraction  of  one  (Yi  =  1)  is  maintained  as  a  nontrivial  input  ‘control’  at  the  inlet  (ri) 
for  two  seconds,  followed  by  a  clearing  pulse  of  one  second  with  no  input  ‘control’  at  the 
inlet  (Fig.  2).  At  the  end  of  the  clearing  pulse  a  normalized  phosphine  mass  fraction  of 
one  (Y2  =  1)  is  maintained  as  a  second  nontrivial  input  ‘control’  at  the  inlet  (IT)  for  two 
seconds,  followed  by  a  clearing  pulse  of  one  second  with  no  input  ‘control’  at  the  inlet. 
Time  integration  is  implemented  using  a  backward  Euler  method  with  fixed  time  steps 
(0.02  s).  Intermediate  solution  values  are  stored  at  each  time  step  to  be  used  later  in 
construction  of  the  POD  basis  elements. 


Construction  of  POD  Modes  The  reactant  transport  simulation  described  above 
provides  a  multivariate  data  set  consisting  of  K  vectors  Xn  =  {x^,  x^2 . . . . ,  x^K},  each 
vector  representing  N  nodal  values  mass  fraction  of  the  nth  species  at  different  times 
during  the  pulsing  sequence  This  original  data  set  Xn  is  transformed  to  a  new  set  of 
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uncorrelated  variables  (POD  modes) 


^  =  nM,  •■■>*}  =  va,  (8) 

where  the  columns  of  <f>n  =  (pn‘2-  ■  ■  ■  >  (PnK }  are  the  eigenvectors  of  the  product  ma¬ 

trix  (Xn' Xn)^  =  A nifinii  ranked,  in  descending  order,  with  respect  to  the  associated 
eigenvalue.  The  prime  superscript  denotes  the  transpose  of  the  matrix.  The  POD  modes 
Zn  are  orthogonal  =  \nAj,  and  the  transformation  of  variables  preserves  the 

data  variability 

K  K  K 

^2(Xn'Xn)kk  =  J2(Zn'Zn)kk  =  ^2  Xnk.  (9) 

k= 1  k= 1  k= 1 


Expansion  of  the  original  data  Xn  in  terms  of  the  most  significant  POD  modes  minimizes 
the  mean  square  error  of  a  reduced  basis  representation  [20] 

M 

(Xn)ij  = 

k=  1 

where  M  <  K.  The  most  significant  POD  modes  are  those  corresponding  to  the 
largest  eigenvalues,  since  the  ratio  of  an  eigenvalue  to  the  summation  of  eigenvalues, 
Xnk/Y^f=i  -Wi  gives  the  percentage  of  the  mean  square  error  unaccounted  for  by  elimi¬ 
nating  the  corresponding  POD  mode  z^k  [20]  in  the  reduced  basis  representation. 


(10) 


2.2  The  Control  Problem 


We  seek  to  control  the  mass  fractions  of  TMI  and  phosphine  at  the  inlet  (rj  of  the 
reactor  in  order  to  obtain  a  desired  flux  qT{t )  of  reactants  at  point  xp  on  the  susceptor 
(r5).  The  general  flux  at  a  point  xp  is  given  by 


Q(t) 


qin(t ) 
Qp(t) 


n  WP  oY 2 
2  ii'  2  dn 


\  Xp 


(11) 


where  Wjn  and  Wp  are  the  molecular  weights  of  indium  (a  component  of  Y\  and  Y3 )  and 
phosphorus  (W),  respectively. 


Under  the  dilute  approximation,  the  steady  state  solutions  for  v.  p,  and  T  can  be  used 
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ill  (5).  The  full  Dirichlet  boundary  control  problem  can  then  be  formulated  as 


dYn 

dt 


-f  v  ■  VYn 
Yn{ 0, x) 
Yi(t,x) 
Y2(t,x) 
Y3(t,x ) 
Yn(t,x ) 
dYn(t,  x) 
Bn 

n 


iv  •  (pDnVYn)  +  XnYi 
Vno(x) 


ux{t) 

U2{t) 

0 

0 

0 

1,2,3 


on  n 
on  n 
on  n 
on  T5 

on  T2  U  T3  U  T4  U  T6 


(12) 


where  Yi.Y2.  and  Y3  refer  to  the  mass  fractions  of  TMI,  phosphine,  and  MMI,  respec¬ 
tively,  and  ui(t),u2(t)  are  the  controls  corresponding  to  TMI  and  phosphine,  Ai(T)  = 
-k0e^s/RT\  A 2(T)  =  0,  and  A 3(T)  =  {W3/W1)k0eYE/RTl 


We  consider  a  finite  time  horizon  problem  of  minimizing  the  cost  function 

J ( u ,  y)  =  [  [u-'Ru  +  (q  -  qr)'Q{q  -  qr)\  dt , 

Jo 

where  the  flux  q(t)  tracks  the  desired  flux  qr{t). 


(13) 


Penalty  Boundary  Formulation.  We  use  a  penalty  boundary  formulation  of  the  time- 
dependent  species  equations  (5)  to  describe  the  transport  of  TMI,  phosphine,  and  MMI 
in  the  reactor  for  implementation  of  the  control  problem 


dYn 

dt 


+  v  ■  VYn 
Yn( 0,  x) 

dYijt,  x) 
Bn 

dY2{t,  x) 
Bn 

dY3(t,  x) 
Bn 

9Yn(t,  x) 
Bn 

dYn(t,  x) 

Bn 


n 


iv  •  (PDnVYn)  +  AnTi 

Vno(x) 


(Yi(t,x)  - 

- Mt )) 

on 

n 

C Y2(t,x )  - 

-  Mt)) 

on 

n 

Y3(t,x ) 

on 

n 

Yn(t,  x) 

on 

T5 

on 

T2  U  T3  U  T4  U  T6 

1,2,3 


where,  as  before,  Y1;Y2,  and  Y3  refer  to  the  mass  fractions  of  TMI,  phosphine,  and  MMI, 
respectively,  and  ui(t),u2(t)  are  the  controls.  Under  sufficient  regularity,  one  can  argue 
that  solutions  of  (14),  in  the  limit  as  e  — >  0,  approximate  solutions  for  the  problem 
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with  conditions  Yn(t,x)  =  un(t )  on  Tl,  and  Yn(t,x )  =  0  at  T5  (see  [5,  12]  for  related 
discussions).  A  value  of  e  =  1  x  10-3  is  used  for  the  results  presented  here. 

Writing  (14)  in  weak  form  with  test  functions  w3 ,  we  obtain 
f  ^  =  —  f  (v  •  VW)  Wj  dQ  —  f  DnVYn  ■  Vwj  dtd  +  f  —iv3DnVYn  •  V pdQ 

+  /  A nYnWjdD,  +  -  /  w,DnYn  ds - /  WjDnunds  , 

in  e  Vn,r5  e  ./ri 

(15) 

where  n  =  1,  2,  3  and  u%  =  0.  Spatially  dependent  values  for  v.  T,  p,  and  Dn  are 
interpolated  from  the  nodal  values  obtained  from  the  CFD  simulations. 

Discretization  Method  I:  finite  elements.  Two  discretization  formulations  are  ap¬ 
plied  to  (15).  The  first  formulation,  a  finite  element  approximation,  produces  the  full 
system  (which  plays  the  role  of  a  reactor  simulator).  In  this  case,  the  species  mass  frac¬ 
tions  (15)  are  approximated  using  standard  finite  element  discretization  with  N  =  453 
quadratic  interpolation  functions  yy  and  N  nodal  coefficients  y^d 

N 

Y"(t,  x)  (16) 

i= 1 


Choosing  the  test  functions  to  be  Wj  =  V’j,  3  =  1,  2, . . . ,  IV,  we  obtain  in  a  standard  way 
the  matrix  equation 


y3N(t )  =  A3Ny3N(t )  +  B3Nu(t )  u(t) 


ui(t) 
u2{t)  ’ 


(17) 


where  A3N  is  an  3 TV  x  3 N  matrix,  B3N  is  a  3 N  x  2  matrix,  and  the  control  u  is  a  control 
vector. 


Discretization  Method  II:  POD  modes.  The  second  discretization  formulation  pro¬ 
duces  the  reduced  basis  model.  In  this  case,  we  first  use  the  POD  modes  {z^k}  to  obtain 
the  POD  basis  elements 

N 

^ nk(x )  =  '52(znk)i'lpi(x),  k  =  1,  2,  ...  ,  K,  (18) 

i= 1 

where  the  functions  ipi(x)  are  the  finite  element  quadratic  interpolation  functions  and 
n  =  1,2,  3.  The  mass  fraction  for  the  nth  species  is  approximated  as  a  linear  combination 
of  the  POD  basis  elements  corresponding  to  the  Mn  most  significant  POD  modes 

Mn 

Y™n(t,x)  =J2ym(t)^ni(x)  (19) 

i= 1 


10 


where,  in  this  case,  Mn  «  K  «  N .  Application  of  this  approximation  to  (15)  (in  this 
case  we  use  POD  test  functions  Wi  =  T™ ,  «  =  1,2,...,  M)  yields 


yM(t)  =  AMyM(t )  +  BMu(t)  u(t ) 


Ui(t) 
U2{t)  ’ 


(20) 


3 

where  M  =  V)  Mn .  AM  is  an  M  x  M  matrix,  and  BM  is  an  M  x  2  matrix. 

n=l 

Tracking  Control.  To  control  the  reduced  order  system  (20),  we  observe  the  flux  q(t) 
of  indium  and  phosphorus,  to  the  center  of  the  substrate  at  xp  (11).  The  fluxes  are 
approximated  as 


qM(t ) 


(- HM)'yM(t ) . 


ggi  k 

dn 


yikit)  +  D3 

Xp 


wIn 

w3 


m3 


E 


V2  k(t) 

Xp 


d^sk 

dn 


V3k(t) 

Xp 


(21) 


We  seek  the  optimal  control  u*  for  (20)  such  that  the  output  qM  tracks  a  signal  qT  (the 
desired  flux  at  xp).  minimizing  the  generalized  performance  index  ( e.g .,  see  [3,  36]) 


u'Ru  +  (qM)'Q2qM  +  (qM 


qTyQi(qM 


qT)  dt  , 


(22) 


where  Q2  keeps  qM  bounded,  qM  =  ( HM)'yM ,  (HM)'  =  I  —  LM(HM )', 

LM  =  HM ((HMy HM)~l ,  and  I  is  the  identity  matrix.  Choosing  Q2  =  r2I ,  R  =  /,  and 


Qi 


Hi  0 
0  r  22 


(23) 


we  can  rewrite  this  performance  index  as 


J(yo,u(-))  =  [ 

Jo 


[u'Ru  +  ( yM  -  y™)'Q(yM  -  y™)\  dt ,  (24) 

where  r n,  r22,  and  r2  are  design  parameters,  Q  =  HMQ2(HM)'  +  HMQ\{HM)'  1  and 

V${t)  =  LqT(t)  (25) 


is  the  desired  state  trajectory.  This  is  a  standard  formulation  for  a  ‘tracking’  control 
problem  [3,  36]  for  which  a  complete  theory  is  known  (for  a  summary  and  references  see 
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Chapter  7  of  [8]).  For  the  pulsed  sources  we  specify  time  dependent  target  fluxes  such 
that  the  output  qM  tracks  the  signal 


Qr(t) 


Q'l'zfz  (t) 


with  constants  qr  1,  qrz  and  target  flux  profiles  f\{t)  and  f2(t)  (Fig.  3). 


The  optimal  control  is  given  by 


u 


* 


u\(t) 
u*2(t )  _ 

-KMyM  -R-\BM)'b(t), 


(26) 


(27) 


where  the  gain  is  given  by  (see  p84-85  of  [3],  [36],  or  [8]) 

KM  =  RT1(BMy  n  (28) 

and  II  satisfies  the  algebraic  Riccati  equation  (ARE) 

0  =  IL4m  +  (AM)'U  -  UBMR-1(BM),U  +  Q.  (29) 

The  tracking  variable  b(t)  is  given  by 

b(t)  =  ~(Am  -  BMKM)'b(t )  +  Qy¥(t),  b(T)  =  0,  (30) 

where  T  =  Tp  +  A,  Tp  is  the  total  elapsed  time  for  the  TMI  (or  phosphine)  pulse  cycle, 
and  A  is  five  times  the  dominant  time  constant  associated  with  the  eigenvalues  of  AM  — 
BmKm  (see  p86  of  [3]).  Since  there  is  a  delay  between  application  of  the  control  (mass 
fraction  of  sources  at  the  inlet)  and  detection  of  the  corresponding  signal  (flux  at  the 
center  of  the  substrate),  the  time  dependent  tracking  variable  plays  a  crucial  role  in  the 
optimal  control,  by  allowing  it  to  anticipate  and  follow  the  time  varying  target  fluxes. 

State  Estimation.  Application  of  the  tracking  control  to  the  reduced  order  model  (20) 
yields 

yM  =  AMyM -BMKMyM -BMR-l(BM)'b(t) 
q(t )  =  (. HM)’yM(t )  1  j 

with  KM  and  b(t)  defined  in  (28)  and  (30),  respectively.  Note  that  this  formulation 
requires  complete  reduced  order  state  feedback  for  the  gain  KM  of  (28).  Since  the  full 
state  of  the  system  (or  even  the  reduced  order  system)  in  the  reactor  can  not  be  observed, 
we  implement  a  state  estimator  or  compensator  design  based  on  observation  (21)  of  the 
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Figure  3:  Indium  and  phosphorus  target  flux  profiles. 


flux  at  the  center  of  the  substrate.  This  yields  (see  [36]  or  Chapter  8  of  [8]  and  the 
references  therein)  the  coupled  system 

yM  =  AMyM  -  BMKMyf  -  BMR~\BM)'b{t) 
yeM  =  Acyf  +  FM(HM)'yM  -  BMR~1(BM)'b(t),  ^ 

where  y ^  is  the  Mxl  vector  approximation  to  the  approximate  state  yM ,  the  compen¬ 
sator  system  operator  Ac  is  given  by 

Ac  =  Am  -  Fm{Hm)'  -  BmKm,  (33) 

and  the  compensator  gain  (Kalman  filter)  is  given  by 

Fm  =  Y,HmV~1.  (34) 

The  matrix  E  satisfies  the  dual  ARE  given  by 

Am E  +  E (AM)'  -  T,HmV~1  (Hm)'H  +  U  =  0,  (35) 

where  U  and  V  are  symmetric  positive  semi-definite  and  symmetric  positive  definite 
matrix  design  parameters,  respectively.  We  choose  U  =  I  and  V  —  r%,  where  r 3  is  a  third 
parameter  at  our  disposal  in  designing  the  overall  feedback  control/compensator  system. 
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This  implementation  of  the  state  estimator  yields  the  following  closed-loop  system  for 
‘optimal’  control  of  the  reduced  order  model 


with  the  optimal  control 

u*  =  - KMy f  -  R-\BM)'b{t).  (37) 

Control  of  the  Full  System.  Use  of  the  design  given  in  (36)  and  (37)  can  be  expected  to 
produce  a  stabilized  and  generally  efficient  system  for  control  of  the  reduced  order  model 
(20).  However,  this  is  not  the  issue  of  practical  importance  in  our  efforts.  The  goal,  of 
course,  is  to  design  a  feedback  control/compensator  system  for  (14)  (or  in  weak  form 
(15)),  which  we  hope  will  be  a  good  approximation  (for  proper  choice  of  e)  for  the  actual 
physical  dynamics  (l)-(5),  i.e.,  transient  solutions  of  (5)  with  (l)-(4)  in  steady  state. 
Thus,  the  real  measure  of  the  value  of  the  control/compensator  system  designed  above 
via  the  reduced  order  model  is  how  well  it  performs  when  used  in  the  actual  physical 
system.  Short  of  applying  the  control/compensator  system  to  the  physical  reactor  in 
experiments,  our  best  assessment  of  its  utility  is  when  applied  to  the  ‘full’  system  (17). 
That  is,  we  should  computationally  test  the  reduced  order  control/compensator  design 
based  on  (28)-(30),  (33)-(37)  in  the  system  (17). 

Application  of  the  reduced  order  tracking  control  to  the  full  system  with  the  reduced 
order  state  estimator  yields 


where  HN  is  the  observation  vector  for  the  full-dimensional  system  qN  =  (HN)'yN(t), 
and  HM ,  KM ,  gM ,  and  FM  are  as  defined  in  (21),  (28),  (30),  and  (34),  respectively.  The 
optimal  control  is  given  by  (37).  We  report  on  our  simulations  for  system  (38)  with 
different  values  of  the  design  parameters  rn  and  r22  hi  the  results  of  Section  3.  Values  of 
rg  and  were  varied  and  then  fixed  at  nominal  performance  values  in  the  calculations 
reported  below. 

Implementation.  The  governing  equations  (14)  are  nondimensionalized  prior  to  the 
calculation  of  the  coefficient  matrices  (Lq  —  lx  10-2  m,  Dq  =  1.15  x  10-3  m2/s,  and 
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p0  =  4.04  x  10-2  kg/rn3).  The  POD  modes  and  the  coefficient  matrices  AN ,  BN ,  HN . 
AM ,  BM ,  and  HM  in  (38)  are  calculated  using  in-house  C  programs.  All  other  matrix 
calculations  are  implemented  in  the  Matlab  computing  environment  (The  Math  Works 
Inc.,  Natick,  MA);  the  optimal  gain  matrix  KM  and  Riccati  solution  II  are  determined 
using  Matlab’s  lqr()  function. 

The  dynamical  equations  are  integrated  using  a  semi-implicit  extrapolation  integration 
method  for  stiff  equations  with  adaptive  step  size  control  ( stifbsQ ,  and  simprQ  in  [35]). 
During  the  clearance  cycles  the  optimal  control  (mass  fractions  at  the  inlet)  may  become 
negative.  Since  negative  mass  fractions  are  not  physically  meaningful,  we  implement  a 
truncation  procedure  at  each  step  of  the  integration,  which  sets  negative  control  values 
equal  to  zero.  The  resulting  control  is  suboptimal,  but  it  will  be  seen  that  control 
authority  is  maintained. 
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3  Results 


We  report  in  sequential  form  results  from  the  series  of  computations  and  simulations 
described  in  detail  in  Section  2. 


3.1  CFD  Simulation  Results 

We  first  report  on  results  of  the  CFD  package  simulations  described  in  Section  2.1. 
Contour  plots  of  the  steady-state  solutions  of  the  temperature  and  the  x-component  of 
the  velocity  (Fig.  4a-b)  give  an  indication  of  the  steady-state  transport  conditions  in 
the  reactor.  With  the  top  wall  and  inlet  maintained  at  room  temperature,  there  is  a 
steep  temperature  gradient  (Fig.  4a)  upstream  from  the  substrate.  In  the  region  above 
the  substrate,  the  isotherms  run  parallel  to  the  hot  substrate  and  the  opposite  cold  wall. 
The  velocity  of  the  gas  increases  more  than  4-fold  as  it  passes  in  the  vicinity  of  the  hot 
substrate  (Fig.  4b),  while  the  y-component  of  the  velocity  (not  shown)  remains  small 
throughout  the  reactor. 


(a)  Temperature 


(b)  X-component  of  velocity 


Figure  4:  Contour  plot  of  steady-state  values  for  the  (a)  temperature  in  100°K  steps  with  Tmoa:=750oK 
and  T,„jn=350°K  ,  (b)  x-component  of  the  velocity  in  0.1  m/s  steps  with  vmoa:=0.45  and  vm.;?)=0.05 
m/s. 
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A  plot  of  the  flux  of  In  and  P  at  the  substrate  center  as  a  function  of  time  (Fig.  5)  shows 
that  the  flux  plateaus  approximately  1  s  after  the  reactant  is  introduced  into  the  reactor 
and  is  not  completely  cleared  from  the  reactor  until  one  second  after  the  cessation  of  the 
reactant  pulse. 


0  5  10  15  20 


NONDIMENSIONAL  TIME 

Figure  5:  Flux  of  In  and  P  to  the  substrate  center  as  a  function  of  time. 
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3.2  Construction  of  POD  Modes 


The  POD  modes  for  each  species  are  constructed  from  150  snapshots  ( K  =  150)  taken 
during  the  three  second  cycle  (2  s  pulsing,  1  s  clearing)  of  each  source  species.  Each  solu¬ 
tion  vector  represents  the  species  mass  fraction  at  the  453  nodal  points  and  corresponds 
to  a  time  increment  of  0.03-s  in  the  time  range  from  0  to  3  seconds.  A  plot  of  the  cap¬ 
tured  variability  i  Ary /  J2k  ^nk)  as  a  function  of  the  number  of  modes  (Mn  <  100) 
used  (Fig.  6)  shows  that  the  original  data  is  well-represented  by  9  modes,  or  fewer.  This 
strongly  suggests  using  Mn  <  9  in  our  reduced  order  model. 


Scree3.m  26-Apr-1999 


Figure  6:  Total  percent  variability  captured  as  a  function  of  the  number  of  modes  for  each  species 
(TMI,  MMI,  phosphine). 


While  the  above  comments  suggest  the  proper  order  for  accurate  reduced  order  model 
simulations,  there  are  additional  order  questions  related  to  the  control  system  to  be  used 
in  determining  reduced  order  gains  and  compensators.  The  ranks  of  the  controllability 
and  observability  matrices  have  sometimes  been  found  to  be  useful  criteria  (see  the  dis¬ 
cussion  in  [9,  10,  26])  to  help  in  the  choice  of  the  number  of  modes  to  use  in  control 
design  applications  for  the  reduced  basis  representation  (19).  The  controllability  of  the 
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linear  system  (20)  is  determined  from  the  rank  of  the  controllability  matrix  C,  where 


C(AM,  Bm)  =  [Bm  |  AmBm  |  (Am)2Bm  |  ...  |  (AM)M-1BM],  (39) 

Using  standard  results  from  optimal  control  theory  [17],  the  M-dimensional  linear  system 
(20)  is  controllable  if,  and  only  if,  the  rank  of  the  controllability  matrix  is  equal  to  M. 
Similarly,  the  rank  of  the  observability  matrix  O  indicates  the  observability  of  the  linear 
system  (31),  where 

0(Am,Hm)  =  [Hm  |  (Am)'Hm  |  ((Am)2),Hm  |  ...  |  ((AM)M_1)'HM].  (40) 

For  the  system  under  investigation  in  this  paper,  ranks  of  the  controllability  and  observ¬ 
ability  matrices  are  recorded  in  Table  1  as  a  function  of  the  number  of  modes  in  the 
reduced  order  model.  Based  on  these  results,  which  suggest  that  the  maximum  ranks  of 
the  controllability  and  observability  matrices  are  ten  and  twelve,  respectively.  However, 
we  should  remember  that  these  ranks  are  for  the  reduced  order  system  and  we  are  ap¬ 
plying  the  control  in  the  full  order  system.  Nonetheless,  we  first  used  ten  POD  modes 
(M  =  10  :  M\  —  4,  M2  =  4,  M3  =  2)  in  our  control/compensator  applications  to  construct 
our  reduced  order  model  (20)  and  obtained  adequate  control  authority.  With  further  ex¬ 
perimentation,  we  found  that  by  adding  more  modes  (M  =  19  :  Mi  =  8,  M2  =  8,  M3  =  3), 
though  no  longer  being  assured  of  the  controllability  and  observability  of  the  reduced  or¬ 
der  system,  we  were  able  to  achieve  tighter  control  and  reduce  the  oscillations  of  the 
system  as  it  tracks  the  desired  flux  values  (see  Fig.  9  below)  in  the  full  order  system. 


Table  1:  Rank  of  the  Controllability  and  Observability  Matrices 


Reduce  Order  Model 

Dimension  M  and  %  Variability  Represented 

TMI  Phosphine  MMI  Rank 


Mi 

% 

m2 

% 

m3 

% 

C 

o 

5 

(99.930) 

5 

(99.942) 

3 

(99.995) 

10 

12 

5 

(99.930) 

5 

(99.942) 

2 

(99.929) 

10 

12 

4 

(99.756) 

4 

(99.793) 

2 

(99.929) 

10 

10 

5 

(99.930) 

4 

(99.793) 

2 

(99.929) 

10 

11 

5 

(99.930) 

3 

(99.264) 

2 

(99.929) 

9 

10 

5 

(99.930) 

4 

(99.793) 

1 

(98.230) 

10 

10 

8 

(99.999) 

8 

(99.999) 

3 

(99.995) 

10 

12 
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3.3  Control  of  Full  System  Using  the  Reduced  Order  State  Es¬ 
timator 

A  set  of  representative  results  are  depicted  here,  where  the  reduced  order  state  estimator 
is  used  to  control  the  full  system  (38)  using  non-dimensional  tracking  values  of  qn  = 
0.0855  and  qr 2  =  0.15.  In  this  example,  we  fix  the  design  parameters  rn  =  5000, 
r2  —  1  x  10-6,  and  r3  —  1  x  104,  and  vary  the  design  parameter  r22  =  150,200,  or 
300.  Initially,  the  solution  vector  for  the  estimated  state  is  given  a  small  nonzero  value 
(y^1  =  1  x  10— 4 ) ,  while  the  full  system  solution  is  given  an  initial  value  of  zero  ( y3N  =  0), 
where  N  =  453  and  M  =  19  (Mi  =  8,  M2  =  8,  M3  =  3).  The  tracking  values  and  pulse 
profiles  were  chosen  so  that  the  integrated  quantity  of  material  deposited  represents 
approximately  one  monolayer  each  of  indium  and  phosphorus  (assuming  a  100  surface 
with  4.1  x  1018  sites  per  m2  [40]).  Because  of  the  truncation  employed  (see  Section  2.2) 
the  results  shown  here  are  not  optimal.  However,  we  demonstrate  that  the  reduced  order 
state  estimator  or  compensator  system  is  still  capable  of  substantial  control  authority 
when  applied  to  the  full  system. 
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The  controls  (ui(t) ,  u2(t))  are  plotted  in  Fig.  7  as  a  function  of  time  along  with  the 
(nondimensional)  target  flux  pulse  cycle  (indicated  by  the  dotted  lines).  From  Fig.  7  it 
can  be  seen  that  the  control  rises  and  falls  in  advance  of  the  tracking  pulses.  Recall  that 
this  is  necessary  because  of  the  delay  between  introduction  of  reactant  at  the  inlet  and 
transport  of  reactant(s)  to  the  substrate.  It  is  also  evident  from  Fig.  7  that  the  control 
values  oscillate  as  the  system  adjusts  to  match  the  observed  flux  to  the  desired  tracking 
flux.  Both  the  peak  values  and  the  magnitude  of  the  oscillations  of  u2  (corresponding  to 
the  phosphine  mass  fraction  at  the  inlet)  increase  as  r2 2  increases,  i.e.,  as  more  weight 
is  placed  on  achieving  the  desired  flux  and  less  weight  is  placed  on  the  cost  of  control. 
As  expected,  the  profile  for  u\  (corresponding  to  the  TMI  mass  fraction  at  the  inlet)  is 
unchanging  since  r n  is  held  constant  in  this  example. 


Figure  7:  Control  values  (TMI  and  phosphine  normalized  mass  fractions  at  the  inlet)  as  a  function 
of  time  for  different  values  of  the  control  parameter  r22-  150  (solid  line),  200  (dashed  line),  and  300 
(dash-dot  line).  A  nondimensional  representation  of  the  target  flux  profile  is  also  shown  (dotted  line) 
for  reference.  The  control  rises  and  falls  in  advance  of  the  tracking  pulses. 
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A  plot  of  the  observed  flux  as  a  function  of  time  (Fig.  8)  shows  that  the  system  is  able 
to  closely  track  the  time  dependence  of  the  desired  flux  profile  without  significant  delays. 
The  ability  of  the  system  to  match  the  target  flux  is  sensitive  to  the  design  parameter 
t 22  with  close  agreement  for  the  case  of  r 22  =  200. 


Figure  8:  Observed  fluxes  as  a  function  of  time  for  different  values  of  the  control  parameter  r22:  150 
(solid  line),  200  (dashed  line),  and  300  (dash-dot  line).  The  target  flux  profile  is  also  shown  (dotted  line) 
for  reference.  The  system  is  able  to  closely  track  the  desired  flux  profile. 


Some  oscillation  of  the  observed  flux  can  be  seen  in  Fig.  8.  The  amplitudes  of  these 
oscillations  increase  with  increasing  control  (larger  values  of  r2 2  or  rn)  and  increase  as 
the  size  of  the  reduced  order  system  is  decreased.  A  comparison  of  the  observed  flux 
using  two  different  reduced  order  dimensions  (Fig.  9)  shows  that  reducing  the  order  from 
M  =  19  to  M  —  10  (the  apparent  rank  of  the  controllability  matrix  for  the  reduced  order 
system)  increases  the  amplitude  of  the  oscillations. 

The  ability  of  the  system  to  track  the  pulse  turn  off  is  limited  when  the  target  pulse  fall 
off  is  steep.  Since  there  is  no  mechanism  for  removing  reactant  from  the  system  other 
than  the  carrier  gas  transport  and  absorption  at  the  substrate,  the  observed  flux  can  not 
fall  off  faster  than  some  characteristic  rate  associated  with  the  physical  parameters  of 
the  problem  ( e.g .,  reactor  dimension  and  carrier  gas  flow  rate).  As  a  result,  the  observed 
flux  may  have  a  tail  that  persists  substantially  beyond  the  target  pulse  drop. 
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Figure  9:  Comparison  of  observed  flux  as  a  function  of  time  for  two  reduced  order  dimensions:  M  =  19 
(dashed  line)  and  M  =  10  (solid  line).  The  target  flux  profile  is  also  shown  (dotted  line)  for  reference. 
The  design  parameter  values  for  these  simulations  were  as  follows:  qxi  =  0.0855,  qx i  =  0.15,  r n  =  1600, 
r 22  =  500,  r-2  —  1  x  10-4,  and  r%  =  1  x  104.  The  magnitude  of  the  oscillations  increases  when  the  order 
is  reduced. 

Fig.  10  plots  the  norm  of  the  difference  between  the  full  system  solution  and  the  state 
estimated  solution  as  a  function  of  time.  The  error  decreases  from  an  initial  non-zero,  the 
value  of  which  depends  on  the  initial  conditions,  then  begins  to  rise  as  the  control  values 
rise  in  anticipation  of  the  rising  edge  of  the  target  flux  pulse.  The  error  decreases  again 
as  the  target  flux  plateaus  and  the  control  values  are  relatively  constant,  then  rises  again 
as  the  control  values  fall  in  anticipation  of  the  target  flux  drop  off.  The  error  reaches 
a  minimum  again  during  the  clearance  time  period,  after  which  the  cycle  of  rising  and 
falling  is  repeated  again  for  the  phosphine  pulse. 
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Figure  10:  The  norm  of  the  difference  between  the  full  state  and  the  POD  estimated  state  for  different 
values  of  the  control  parameter  150  (solid  line),  200  (dashed  line),  and  300  (dash-dot  line).  A 
nondimensional  representation  of  the  target  flux  profile  is  also  shown  (dotted  line)  for  reference.  The 
error  rises  as  the  control  values  change  and  decreases  for  during  periods  of  (relatively)  constant  control. 

4  Conclusion 


We  have  demonstrated  a  computational  implementation  of  reduced  order  feedback  control 
of  pulsed  HPCVD  III-V  film  growth  involving  the  transport  of  multiple  species  with 
linear  gas  phase  reactions.  We  implement  feedback  control  using  a  reduced  order  state 
estimator  based  on  observations  of  the  fluxes  of  the  Group  III  and  Group  V  reactants 
at  the  substrate  center.  These  observations  are  compatible  with  current  PRS  sensing 
technology.  The  controls  are  chosen  so  that  the  output  fluxes  track  time  dependent 
target  fluxes,  similar  to  the  pulsed  sources  currently  employed  for  HPCVD  film  growth. 
Because  the  control  values  must  be  constrained  to  be  positive  (positive  mass  fraction) 
the  resulting  truncated  control  is  suboptimal,  but  the  reduced  order  model  design  is  still 
capable  of  substantial  control  authority. 

Use  of  the  POD-based  design  method  allows  us  to  reduce  the  order  of  the  system  with 
respect  to  a  standard  finite  element  representation,  from  3 N  =  2159  to  M  =  19.  These 
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results  suggest  that  real-time  feedback  control  with  partial  state  observations  is  a  feasi¬ 
ble  goal  for  HPCVD  reactors  operating  in  steady  state  flow  regimes  with  pulsed  vapor 
reactant  inputs.  The  positive  results  on  this  linear  system  suggest  one  possible  approach 
for  treating  more  complex  situations  that  may  be  encountered  when  nonlinear  gas  phase 
or  surface  phase  reactions  are  present:  linearizing  the  equations  about  an  optimal  open 
loop  solution  and  then  using  a  reduced  order  feedback  design  on  the  linearized  system. 
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